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By using the point model of reaction kinetics we have studied the stochastic properties of the 
lifetime of small systems controlled by autocatalytic reaction A + X ^ X + X, X ^ B. Assuming 
that a system is living only when the number of autocatalytic particles is larger than zero but 
smaller than a positive integer A'^, we have calculated the probability of the lifetime provided that 
the number of substrate particles A is kept constant by a suitable reservoir, and the end-products B 
do not take part in the reaction. We have shown that the density function of the lifetime is strongly 
asymmetric and in certain cases it has a well-defined minimum at the beginning of the process. It 
has been also proven that the extinction probability of systems of this type is exactly 1. 
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I. INTRODUCTION 

Finding models for calculations of lifetime probability of systems evolving in accordance to well-defined rules seems 
to be a very attractive problem 0, 0, S S 13 not only in physics but also in many other fields of sciences. The 
system evolution is usually controlled by stochastic transformations of components which are called in the sequel 
" particles" . In the present model the transformation is realized by autocatalytic reactions which play important role 
in organization and coordination of systems having small spatial volume and containing few particles only. 

In this paper we are dealing with small systems which are in living state when the number of autocatalytic particles 
remains within the interval (0, A'^), where A'^ is a positive integer, while the number of substrate particles is kept 
constant by a suitable reservoir. One can say that both the absence and the large "dose" of X particles lead to the 
extinction of the system. We would like to determine the lifetime probability of such systems. 

We choose from many possible autocatalytic reactions one of the simplest, namely A + X ^ X + X ^ B, where 
A denotes the substrate particles the number of which is assumed to be constant, and X symbolizes the autocatalytic 
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FIG. 1: Illustration of a reversible autocatalytic reaction by a simple graph. 

particles. The end-product particles B do not take part in the reaction. For the sake of descriptiveness the reaction 
A+X ^ X+X is illustrated by a simple graph shown in FIC^ where the white circle corresponds to an autocatalytic 
while the black one to a substrate particle. 




FIG. 2: Illustration of a system transformation caused by autocatalytic reaction. 

A system consisting of A'^ substrate and Nx autocatalytic particles can be represented by a set of complete graphs 
if AT on = Na + Nx points, where Na is fixed and max A^x is denoted by A^. A randomly chosen point of a graph 
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may be either A or X, the only requirement is that the total numbers of A and X points must be Na and Nx, 
respectively. Clearly, the number of complete graphs characterizing a given systems is given by 1/B{Na + 1, Nx + 1), 
where B{x, y) is the beta function. Since every two points adjacent, the possibility of the interaction between any two 
particles of the system is the same. Therefore, it is obvious that the point model of reaction kinetics introduced by 
[g| can be applied. As an example, FIG. |2] illustrates a transformation which brings about a system oiV — Q particles 
from a system of = 5 particles and vice versa. 

In Section II we derive equations determining the probability n) of finding n autocatalytic particles at the time 
moment t in a system of volume V provided that at t = the number of particles was i. In Section III we show that 
the extinction probability of the system is 1, while in Section IV results of exact calculations for = 3 and N — h 
are presented. The main conclusions are summarized in Section V. 



II. SYSTEMS WITH LIMITED NUMBER OF AUTOCATALYTIC PARTICLES 

Now, we formulate the basic equation. Let ^(t) be the number of autocatalytic particles X at time instant t > 
in a system open for substrate particles A. In the sequel, we would like to define a special system existing only when 
^(t) remains within the interval (0, N), where iV is a positive integer. We say the system is in the state iS„, if ^(t) = n. 
Obviously, a system will be annihilated when it enters into the state either iSo or Sn- For the sake of simplicity we 
call a system of this type SL system. 

In order to characterize the stochastic behavior of an SL system we define the transition probabilities ^ 

^UW = n\m ^i}^ p,,„it) = p„(t), n = 0, 1, 2, . . . , iV, 
which are determined by the following equations: 

dpo{t) 

— ^=7PiW, (1) 



Pn{t) 
dt 



= -n{NAa + nP - P + j)p„{t) + (n - l)NAapn-i{t) + (n + l){n/3 + j)p„+i{t), 



n = l,2,...,iV-2. 

Since there is no way out from the state Sn the equation for pN-i{t) should have the form 
PJV-i(i) 



dt 



and, finally the equation 



[N - l)iNAa + Nf3-2f3 + 7W-i(t) + (N ~ 2)NAapN-2{t), 

PN{t) 



dt 



{N -l)NAapN-iit) 



determines the probability of the transition Si — » Sn during the time t > 0. Introducing the vector 

pit,N^l) 



( PN-lit) \ 
P2N^2{t) 



V Plit) J 



and the matrix 



/~Dn^i Bn-1 

Cn-2 —Dn-2 Bn-2 









V 



Cn- 






-Dn-3 Bn-3 









2^ 





C2 — -D2 B2 
Ci -dJ 



(2) 



(3) 



(4) 



(5) 



^ In order to simplify the notations the index referring to the initial state Si will be omitted where it does not cause confusions. 
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we can write Eqs. ^ and Q in the following compact form: 

^^H(!:^.;,A._,p(.,A.-l), (6) 

with initial conditions Pn(0) — Sn.i for 1 < n < — 1. The elements of the matrix Aat^i are given by 

Bn — {n ~ l)a, Dn = n{a + b + n — 1), and C„ = (n + + n), (7) 



1< 72 < N - 1, 



where 



a^NA^, and b = ^. (8) 
P P 

One has to note that for the complete description of the process we need also the equations ^ and The matrix 
An-i is a normal Jacobi matrix. (See Ref. '^^^ be proven that the eigenvalues o/ Ajv-i are different, 

negative real numbers. The proof is given in Appendix A. Consequently we obtain 



but 



limp„(t) = 0, if n= 1,2,...,A-1, 



lim po{t) — wq, lim pAr(t) = wn and wq + Wn = 1, (9) 

t — >oo t — ^oo 



i.e. if t — > oo, then an SL system can be found either in Sq or in Sn state. 

III. LIFETIME PROBABILITY 

The random time 6i, 1 < i < N — 1 due to the transitions either Si Sq or St Sn is called the lifetime of the 
SL system. In other words, the lifetime is a random time interval at the end of which the number of autocatalytic 
particles becomes either zero or N. Formulating more precisely we can write 

{9, <t} = {m = o\m = u {m = N\ao) = 

and since {^{t) = 0|C(0) = i} Ci {^{t) = N\^{0) = = 0, we have 

V{e, <t}^ H,{t) = p,,oW + P^,N{t). (10) 
By using Eqs. lO we see immediately that 

\miH,{t) = l, Vl<i<A^-l, (11) 

t — >oo 

i.e. the probability that the lifetime of the system is finite is equal to 1. ^ 

IV. SYSTEMS WITH SMALL NUMBER OF STATES 

The mathematical properties of the equation ^ are well-known Here, we do not wish to make detailed analysis 
of Eq. ||HJ) and do not derive its formal solution, but instead we try to study the basic properties of exact solutions 
in those cases when the maximal number of autocatalytic particles is A'^ = 3 and N — 5, i.e. the number of states is 
four and six, respectively. 



^ This statement is equivalent to that the probability of the extinction of the system is equal to 1. 
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A. Four states systems 

In order to have an insight into the dynamics of the process producing random transitions between different states 
of the system we are deaUng now with a very simple system having only four states, namely Sq, Si, S2, and S^. In 
this case the equations to be solved are 



dpojt) 

dt 
dpijt) 

dt 
dp2{t) 

dt 
dpsjt) 

dt 



IPiit), (12) 

~{aNA + l)pi{t) + 2{p + ^)p2{t), (13) 

aNAPi{t)-2{aNA + l3 + i)p2(t), (14) 

2aNAP2{t). (15) 



The initial conditions to be used are po(0) = P2(0) = P3(0) — and pi(0) = 1. Introducing the notations ((HJ, after 
elementary calculations one obtains 

, d — a — b — 2 , , „ ^ ^ ,c! + a + & + 2, , , ^ 

Po{t) = b (1 - exp{-a;i/30) + b (1 - exp{-c^2/3<}), (16) 

2d uj\ Ad L02 

,.\ cxp{-tj2/3i} + exp{-u;i/3t} a + 6 + 2 
pi{t) = + — (exp{-w2/30 - exp{-wi/3i}) , (17) 

P2{t) = ^ (exp{-co2/30 - exp{-wi/3i}) , (18) 



and 



where 



/ 1 - exp{-W2/3i} 1 - exp{-tJi/3i} ' 
P3(i) = 2 — , (19) 

d \ UJ2 ^1 



1 



while 



d = i/a2 + (6 + 2)2 + 2a(5& + 6), wi = -(A + d) and ^2=2^^"'^)' (^^^ 

A ==2 + 3 (a + 6). (21) 

According to the theorem on eigenvalues nonnegative eigenvalue cannot be occurred, therefore \ > d. Naturally, 
this inequality can be simply proven for all positive a and b by using Eqs. H20|l and (|21|l . ^ 

The probability density functions of the system lifetime 61 are plotted on FIG. 13 at different values of the number 
of A particles (part B of the figure) and at different 7 / (3 ratios (part A of the figure) . It is remarkable that the density 
functions are strongly asymmetric, and the most probable lifetimes measured in units are found near the zero of 
the time parameter. 

As seen from Eqs. H12|l - H15f) we obtain the following limit probabilities: 

, , b (1 + a + b) , , 

t^co a^ + o(l + a + o) t^oo a^ + b{\ + a + b) 

and 

lim pi{t) = lim P2{t) — 0. 

t — >oo t — *oo 

According to this, after elapsing sufficiently long (i.e. infinite) time, the system can be found almost surly in "dead" 
state cither with probability wq in the state So or with probability in the state 53. Obviously wo + W3 = 1. 



Suppose that \ = d, i.e. 

3(a + bf - + (6 + 2f + 2a{5fe + 6) + 2 = 0. 
It is elementary to show the roots of this equation for b to be complex, if < a < 0.4467, and negative, if a > 0.4467. 
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FIG. 3: Probability density functions of the system lifetime 8i. Part A shows the density functions at different 7//? values when Na = 5, 
while part B illustrates the effect of the variation of the number of A particles at fixed 7//? value. The ratio a/P = 1. 



B. Six states systems 

Now, we would like to see how the number of possible living states effects on the properties of the lifetime. If the 
number of states is six, i.e. = 5, then the system has four living and two dead states. According to the notation 
introduced earlier, the living states are iSi, 1S2, iSs, 54, while the dead ones iSo and 55. 
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FIG. 4: Dependence of probabilities pi{t),P2{t),P3{t) and the system lifetime distribution H\(t) on t measured in /3 ^ units. The ratio 
a/P = 1. 



By using numerical method for solving the equation system C3)-(0J in the case of = 5 we obtain the probabilities 
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FIG. 5: Density function of the system lifetime at different values of the number of A particles at 7//3 = 1. The number of living states 
is equal to 4. The time t is measured in /3~^ units. The ratio a/13 = 1. 



Pnit), n = 0, 1, . . . , 5 VS. time parameter curves shown some of them in FIG. 01 What remarkable is the character 
of curves does not change with increasing number of possible states. It can be proven that the higher the number of 
possible states the larger is the probability that the system lifetime is longer than t. 
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FIG. 6: Density function of the system lifetime at different values of the ratio j/lB and at the fixed number of A particles. The number 
of living states is 4 and the time t is measured in l3~^ units. The ratio a/f3 = I. 
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FIG. 7: Density function of the system lifetime at the beginning of the process for three different values of the number of A particles. 
The number of living states is 4, the time t is measured in units and 7//? = 1. The ratio o//3 = 1. 

Density functions of the system lifetime G\ measured in /3~^ units are seen in FIG.|31at values Na = 5, 10, 15. The 
calculations are performed at fixed 7//? = 1 ratio. One can see the curves are strongly asymmetric and the most 
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probable lifetime decreases with increasing Na- The reason for this is clear: the speed of the reaction increases with 
the number of A particles, and so if is large, then the state iSjv can be reached in a time interval shorter than that 
is needed when Na is small. It is remarkable that we see in FIG. El showing the density function of 9i at different 
7//3 ratios. If the decay intensity 7 is approximately equal to the intensity fi of the reverse reaction X + X ^ X + A, 
then a well-defined minimum appears in the density vs. lifetime curve at the very beginning of the process, because 
the autocatalytic production of X particles cannot compensate the effect of the decay process. After elapsing some 
time the role of the autocatalytic process becomes perceptible and, if the ratio 7//? is smaller than 3, then the density 
function reaches a maximum which is followed by a monotonously decreasing part of the curve. The density function 
at the beginning of the process is shown in FIG. for three Na values. 



V. CONCLUSIONS 



Stochastic dynamics of small systema controlled by autocatalytic reaction X + X ^ X + A, X ^ B is studied. 
The dead state of systems is defined and the probability distribution of the lifetime are calculated by exact means. 
The density function of the lifetime is found to be strongly asymmetric (skewed right). The most probable lifetime 
(mode of the lifetime) decreases with increasing number of the substrate particles A. If the self-decay intensity of the 
autocatalytic particles is approximately equal to the intensity of the reverse reaction X + X ^ X + A, then a well- 
defined minimum appears in the density function of the lifetime near the beginning of the process. It is remarkable 
that the basic character of the system's dynamics does not vary significantly with increasing number of possible state. 



APPENDIX A: EIGENVALUES 



We would like to prove that the eigenvalues of the matrix A^-i are different, negative real numbers. Let us 
introduce the matrix 



/ G„ F„ 

Fn-l Gn-1 H„-l 
Gn-2 Hn-2 





V 



\ 





F2 G2 H2 
Fi gJ 



where n e Z^, 



(Al) 



while 



F, = -{j + l)ib + j), Gj=jia + b + j-l), H, = -{j-l)a, 



(A2) 



where a and b are positive real numbers given by Eq. 
which transforms J„ into a symmetric matrix 



Since Fj/Hj^i > there is a diagonal matrix S„ such 



(A3) 



the eigenvalues of which are exactly the same as those of J„. The entries of S„ are 

Si,j = Si J Sj, where Sj 



s,„i, and si — 1. 



It is easy to show that 

/ Gn 







J-1 



Gn — 2 \/ Fn — 'i Hn — 2 



\ 



\ 





y/F^h G2 VWTh. 
VF\H2 Gi j 



(A4) 



and clearly R„ is a symmetric Jacobi matrix. As known, if it is positive definite, then all of its eigenvalues are different, 
real and positive. A symmetric matrix is positive definite if the principal minors of the corresponding determinant 
are all positive. Let us denote by 



=det |R,| 



(A5) 
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the determinant due to the matrix Rj. Taking into account that 2?o = 1 by definition, we have to show that 
Vj > 0, V J = 1, 2, ... , 71. After simple calculation we obtain 

V,=G,V,^^-H,F,^,V,,2, (A6) 

and by induction we have 

which clearly proves the statement. Finally, since the eigenvalues of Rat-i are real, positive and different, it follows 
that the eigenvalues of K^-i are real, negative and different. Q.E.D. 
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